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Abstract. In this paper, we investigate the complex dynamics of a spatial plankton-fish system 
with Holling type III functional responses. We have carried out the analytical study for both one 
and two dimensional system in details and found out a condition for diffusive instability of a locally 
stable equilibrium. Furthermore, we present a theoretical analysis of processes of pattern formation 
that involves organism distribution and their interaction of spatially distributed population with 
local diffusion. The results of numerical simulations reveal that, on increasing the value of the 
fish predation rates, the sequences spots — > spot- stripe mixtures— stripes— > hole- stripe mixtures 
holes— t- wave pattern is observed. Our study shows that the spatially extended model system has 
not only more complex dynamic patterns in the space, but also has spiral waves. 
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1. Introduction 

Complexity is a tremendous challenge in the field of marine ecology. It impacts the develop- 
ment of theory, the conduct of field studies, and the practical application of ecological knowledge. 
It is encountered at all scales (Loehle, 2004). There is growing concern over the excessive and 
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unsustainable exploitation of marine resources due to overfishing, climate change, petroleum de- 
velopment, long-range pollution, radioactive contamination, aquaculture etc. Ecologists, marine 
biologists and now economists are searching for possible solutions to the problem (Kar and Mat- 
suda, 2007). A qualitative and quantitative study of the interaction of different species is important 
for the management of fisheries. From the study of several fisheries models, it was observed that 
the type of functional response greatly affect the model predictions (Steele and Henderson, 1992; 
Gao et al., 2000; Fu et al., 2001; Upadhyay et al., 2009). An inappropriate selection of functional 
response alters quantitative predications of the model results in wrong conclusions. Mathemati- 
cal models of plankton dynamics are sensitive to the choice of the type of predator's functional 
response, i.e., how the rate of intake of food varies with the food density. Based on real data ob- 
tained from expeditions in the Barents Sea in 2003-2005, show that the rate of average intake of 
algae by Calanus glacialis exhibits a Holling type III functional response, instead of responses of 
Holling types I and II found previously in the laboratory experiments. Thus the type of feeding 
response for a zooplankton obtained from laboratory analysis can be too simplistic and misleading 
(Morozov et al., 2008). Holling type III functional response is used when one wishes to stabilize 
the system at low algal density (Truscott and Brindley, 1994; Scheffer and De Boer, 1996; Bazykin, 
1998; Hammer and Pitchford, 2005). 

Modeling of phytoplankton-zooplankton interaction takes into account zooplankton grazing 
with saturating functional response to phytoplankton abundance called Michaelis-Menten models 
of enzyme kinetics (Michaelis and Menten, 1913). These models can explain the phytoplankton 
and zooplankton oscillations and monotonous relaxation to one of the possible multiple equilibria 
(Steele and Henderson, 1981, 1992; Scheffer, 1998; Malchow, 1993; Pascual, 1993; Truscott and 
Brindley, 1994). The problems of spatial and spatiotemporal pattern formation in plankton include 
regular and irregular oscillations, propagating fronts, spiral waves, pulses and stationary spiral pat- 
terns. Some significant contributions are dynamical stabilization of unstable equilibria (Petrovskii 
and Malchow, 2000; Malchow and Petrovskii, 2002) and chaotic oscillations behind propagating 
diffusive fronts in prey-predator models with finite or slightly inhomogeneous distributions (Sher- 
ratt et al., 1995, 1997; Petrovskii and Malchow, 2001). 

Conceptual prey-predator models have often and successfully been used to model phytoplankton- 
zooplankton interactions and to elucidate mechanisms of spatiotemporal pattern formation like 
patchiness and blooming (Segel and Jackson, 1972; Steele and Hunderson, 1981; Pascual, 1993; 
Malchow, 1993). The density of plankton population changes not only in time but also in space. 
The highly inhomogeneous spatial distribution of plankton in the natural aquatic system called 
"Plankton patchiness" has been observed in many field observations (Fasham, 1978; Steele, 1978; 
Mackas and Boyd, 1979; Greene et al., 1992; Abbott, 1993). Patchiness is affected by many fac- 
tors like temperature, nutrients and turbulence, which depend on the spatial scale. Generally the 
growth, competition, grazing and propagation of plankton population can be modeled by partial 
differential equations of reaction-diffusion type. The fundamental importance of spatial and spa- 
tiotemporal pattern formation in biology is self-evident. A wide variety of spatial and temporal 
patterns inhibiting dispersal are present in many real ecological systems. The spatiotemporal self- 
organization in prey-predator communities modelled by reaction-diffusion equations have always 
been an area of interest for researchers. Turing spatial patterns have been observed in computer 
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simulations of interaction-diffusion system by many authors (Malchow, 1996, 2000; Xiao et al., 
2006; Brentnall et al., 2003; Grieco et al., 2005; Medvinsky et al., 2001; 2005; Chen and Wang, 
2008; Liu et al., 2008). Upadhyay et al. (2009, 2010) investigated the wave phenomena and non- 
linear non-equilibrium pattern formation in a spatial plankton-fish system with Holling type II and 
IV functional responses. 

In this paper, we have tried to find out the analytical solution to plankton pattern formation 
which is in fact necessity for understanding the complex dynamics of a spatial plankton system 
with Holling type III functional response. And the analytical finding is supported by extensive 
numerical simulation. 



2. The model system 

We consider a reaction-diffusion model for phytoplankton-zooplankton-fish system where at any 
point (x, y) and time t, the phytoplankton P(x, y, t) and zooplankton H(x, y, t) populations. The 
phytoplankton population P(x, y, t) is predated by the zooplankton population H (x, y, t) which is 
predated by fish. The per capita predation rate is described by Holling type III functional response. 
The system incorporating effects of fish predation satisfy the following: 

§ = rP - - ffgg + d^P, 

(2 I) 

dH _ r> tt /~i H 2 FH 2 i j vy2 tt V ' ' 

with the non-zero initial conditions 

P{x, y, 0) > 0, H(x, y, 0) > 0, {x, y) E Q = [0, R] x [0, R], (2.2) 
and the zero-flux boundary conditions 

OP dH 

— = — = 0, (x, y) e dQ for all t. (2.3) 

on on 



where, n is the outward unit normal vector of the boundary dVl which we will assume is smooth. 

d 2 , <3 2 
dx 2 dy 2 



And V 2 is the Laplacian operator, in two-dimensional space, V 2 = -§s + jpz, while in one- 



dimensional case, V 



2 = d 2 
dx' 2 ' 

The parameters r, Bi, B 2 , D, C\, D\, di, d 2 in model (12.11) are positive constants. We ex- 
plain the meaning of each variable and constant, r is the prey's intrinsic growth rate in the ab- 
sence of predation, Bi is the intensity of competition among individuals of phytoplankton, B 2 is 
the rate at which phytoplankton is grazed by zooplankton and it follows Holling type-Ill func- 
tional response, C\ is the predator's intrinsic rate of population growth, C 2 indicates the number 
of prey necessary to support and replace each individual predator. The rate equation for the zoo- 
plankton population is the logistic growth with carrying capacity proportional to phytoplankton 
density, P/C 2 . D, D x is the half- saturation constants for phytoplankton and zooplankton density 
respectively, F is the maximum value of the total loss of zooplankton due to fish predation, which 
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also follows Holling type-Ill functional response. di,d 2 diffusion coefficients of phytoplankton 
and zooplankton respectively. The units of the parameters are as follows. Time t and length 
x, y 6 [0, R] are measured in days [d] and meters [m] respectively, r, P, H, D and D\ are usually 
measured in mg of dry weight per litre [mg ■ dw/t]; the dimension of B 2 , and C 2 is [a! -1 ], are 
measured in [(mg ■ dw/l^d' 1 ] and respectively. The diffusion coefficients d\ and d 2 are 

measured in [m 2 d~ 1 ]. F is measured in [(mg-dw/Z)^ 1 ]. The significance of the terms on the right 
hand side of Eq.(l) is explained as follows: The first term represents the density-dependent growth 
of prey in the absence predators. The third term on the right of the prey equation, B 2 P 2 / (P 2 + D 2 ), 
represents the action of the zooplankton. We assume that the predation rate follows a sigmoid (Type 
III) functional response. This assumption is suitable for plankton community where spatial mixing 
occur due to turbulence (Okobo, 1980). The dimensions and other values of the parameters are 
chosen from literatures (Malchow, 1996; Medvinsky et al., 2001, 2002 Murray, 1989), which are 
well established for a long time to explain the phytoplankton-zooplankton dynamics. Notably, the 
system (1) is a modified (Holling III instead of II) and extended (fish predation) version of a model 
proposed by May (1973). 

Holling type III functional response have often been used to demonstrate cyclic collapses which 
form is an obvious choice for representing the behavior of predator hunting (Real, 1977; Ludwig 
et al., 1978). This response function is sigmoid, rising slowly when prey are rare, accelerating 
when they become more abundant, and finally reaching a saturated upper limit. Type III functional 
response also levels off at some prey density. Keeping the above mentioned properties in mind, we 
have considered the zooplankton grazing rate on phytoplankton and the zooplankton predation by 
fish follows a sigmoidal functional response of Holling type III. 



3. Stability analysis of the non-spatial model system 

In this section, we restrict ourselves to the stability analysis of the model system in the absence of 
diffusion in which only the interaction part of the model system are taken into account. We find 
the non-negative equilibrium states of the model system and discuss their stability properties with 
respect to variation of several parameters. 



3.1. Local stability analysis 

We analyze model system (12.11) without diffusion. In such case, the model system reduces to 

dP _ p _ o p2 _ B 2 P 2 H 

9H _ n tt r~i H 2 FH 2 ' ' 

~dt - UlU ~~ ° 2 ~F ~ W+Dl 

withP(x,0) > 0,H(x,0) > 0. 

The stationary dynamics of system (13.11 ) can be analyzed from dP/dt = 0, dH/dt = 0. Then, 
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we have 

rP - B 1 P* - ffgf = 0, 

r u r h 2 FH2 - n <y3 ' 2 " > 
Utll - o 2 — - H2+D 2 - u. 

An application of linear stability analysis determines the stability of the two stationary points, 
namely #1(^,0) and E*(P*,H*). 

Notably, E*(P*, H*) presents the kind of nontrivial coexistence equilibria, not only mean one 
equilibrium. In fact, it may be noted that P* and H* satisfies the following equations: 

(r-ff 2 P*)(P* 2 + £ 2 ) _ n FH* H* 
H B~P* -°' H^ + Dl +C2 P^~ Cl -°- 

After solving P* from the second equation above and substituting it into the first equation, we 
can obtain a polynomial equation about H* with degree nine; therefore the equation at least has 
one real root. That is to say, system (13.11 ) may be having one or more positive equilibria, namely 

E*(P*,H*). 

From Eqs. (13.11 ), it is easy to show that the equilibrium point E\ is a saddle point with stable 
manifold in P— direction and unstable manifold in //—direction. The stationary point E*(P*, H*) 
is stable or unstable depend on the coexistence of phytoplankton and zooplankton. In the follow- 
ing theorem we are able to find necessary and sufficient conditions for the positive equilibrium 
E*(P*, H*) to be locally asymptotically stable. 

The Jacobian matrix of model (Q > at E* (P*, H*) is 

/ P*(-B 1 + B 2 H*f 1 (P*,D)) -iltSr 

^ ^ #*(-^ + F/ 2 (iF,L>i)) 

where, 

p*2 _ D 2 H *2 __ D 2 

t I zj* n\ 1 



fi{P m ,D)= , 1M , m , a , f 2 (H*,D) 



Following the Routh-Hurwitz criteria, we get, 

A = P*(B l - B 2 H*h(P*, D)) +H*(^- Ff 2 (H\ D 1 )), 

B = P*H*{B 1 - B 2 H*h{P\ D)) (§£ - F f 2 (H\ D 1 )) + ffi^ 

Summarizing the above discussions, we can obtain the following theorem. 



(3.3) 



Theorem 1. The positive equilibrium E*(P*, H*) is locally asymptotically stable in the PH— plane 
if and only if the following inequalities hold: 

A > and B > 0. 
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Remark 2. Let the following inequalities hold 

B 1 > B 2 H*f x {P%D), C 2 > FP*f 2 (H\D 1 ). (3.4) 

Then A > and B > 0. It shows that if inequalities in Eq. (13.41 ) hold, then E*(P*, H*) is locally 
asymptotically stable in the PH— plane. 

Fig. 1 shows the dynamics of system (13.11) that for increasing fish predation rates at F = (limit 
cycle — no fish), F = 0.02, 0.03, 0.04 (limit cycle) and F = 0.044 (phytoplankton dominance). 
We observe the limit cycle behavior for the fish predation rate F in the range < F < 0.043, and 
for higher values phytoplankton dominance is reported for a non-spatial system (13.11) . 




Phyloplanklon Phyloplankstm 



(a) (b) 

Figure 1: (a) predator-prey zero-isoclines (b) prey zero-isoclines for different value of F with 
parameters value r = 1, B x = 0.2, B 2 = 0.91, d = 0.22, C 2 = 0.2, D 2 = 0.3, D x = 0.1. 



3.2. Global stability analysis 

In order to study the global behavior of the positive equilibrium E*(P* , H*),we need the following 
lemma which establishes a region for attraction for model system (13.11) . 

Lemma 1. The set K = {(P, H) : < P < ^, < H < ^} is a region of attraction for all 
solutions initiating in the interior of the positive quadrant R. 

Theorem 3. If 

P*(r - B 1 P*) 2 < ABirD 2 , rC x E* < D\C 2 B X (3.5) 

hold, then E* is globally asymptotically stable with respect to all solutions in the interior of the 
positive quadrant R. 
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Proof. Let us choose a Lyapunov function 

t 

fv(0 ~ Jh* v 



P t p* fH TT* 

V(P, H) = I q t% + w drj, (3.6) 



where <p(P) = 

Differentiating both sides with respect to t, we get 



dV P-P* dP H — H* dH 

~dt { ' ] ~ P<p(P) ~dt + w H dT' 

Substituting the expressions of ^£ and ^ from Eq. (13.11 ) and putting w = ^^r, we obtained: 

dV _ P-P* ( rP-ffiP 2 rptA / rr U*\2_P*_ .,,( tt u*\2 F{D 2 -H*H) 
~dt ~ —p— \~^(P) 11 ) 1FP~ W \ n ~ U ) (H*+D'(){H**+D'() 

= iJ ^ff- (rP*P - rD 2 - B 1 PP*(P* + P)) - (H - H*) 2 £p 

-,n(TT- U*\2 F(D 2 -H* H) 

w\n n ) { ^ H 2 +D -i^ H ,2 +D 2 i) 
then ^ < for all P, H E K if the condition (EJ) hold. 



4. Stability analysis of the spatial model system 

In this section, we study the effect of diffusion on the model system about the interior equilibrium 
point. Complex marine ecosystems exhibit patterns that are bound to each other yet observed over 
different spatial and time scales (Grimm et al., 2005). Turing instability can occur for the model 
system because the equation for predator is nonlinear with respect to predator population, H and 
unequal values of diffusive constants. To study the effect of diffusion on the system (13.11) . we 
derive conditions for stability analysis in one and two dimensional cases. 



4.1. One dimensional case 

The model system (12.11) in the presence of one-dimensional diffusion has the following form: 

dP _ „ p _ d p2_ B 2 P 2 H , j d 2 P 
dt - ' r D ^ r P 2 +D 2 "1 8x 2 ' 

OH _ r> tt f~t H 2 FH 2 , j d 2 H ^ ' ' 

— _ Ol-tt - C 2 — - H 2 +D 2 i + "2-^2-- 

To study the effect of diffusion on the model system, we have considered the linearized form 
of system about E*(P*, H*) as follows: 

§ = bnU + b 12 V + d l ^, 
f = b 21 U + b 22 V + d 2 ^ 
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where P = P* + U, H = H* + V and 

b n = -P*(B 1 - B 2 H*h(P\ D)), 612 = -B 2 P*/(P* 2 + -D 2 ), 

b 21 = C 2 H* 2 /P* 2 , b 22 = H*(Ff 2 (H*, D l ) - C 2 /P*). 

It may be noted that (U, V) are small perturbations of (P, H) about the equilibrium point 

E*(P*,H*). 

In this case, we look for eigenfunctions of the form 



I n J exp (At + ikx), 



and thus solutions of system (14.21) of the form 



' 71=0 n ' 



where A and k are the frequency and wave-number respectively. The characteristic equation of the 
linearized system (14.21 ) is given by 

A 2 + Pl X + p 2 = 0, (4.3) 

where 

pi = A + (di + d 2 )A; 2 , (4.4) 



p 2 = 5 + d lC / 2 £ 4 + rf 2 P*(5i - B 2 H*h(P\ D)) + dxH\C 2 /P* - Ff 2 (H\ D x )) 



k 2 , (4.5) 



where A and P are defined in Eq. (13 . 3b ■ From Eqs. (14.31 ) — (14.51 ), and using Routh-Hurwitz criteria, 
we can know that the positive equilibrium E* is locally asymptotically stable in the presence of 
diffusion if and only if 

pi > and p 2 > 0. (4.6) 
Summarizing the above discussions, we can get the following theorem immediately. 



Theorem 4. (i) If the inequalities in Eq. (13.41) are satisfied, then the positive equilibrium point E* 
is locally asymptotically stable in the presence as well as absence of diffusion. 

(ii) Suppose that inequalities in Eq. (13.41) are not satisfied, i.e., either A or B is negative or both 
A and B are negative. Then for strictly positive wave-number k > 0, i.e., spatially inhomogeneous 
perturbations, by increasing d\ and d 2 to sufficiently large values, pi and p 2 can be made positive 
and hence E* can be made locally asymptotically stable. 

Diffusive instability sets in when at least one of the conditions in Eq. (14.61) is violated subject 
to the conditions in Theorem 1 . But it is evident that the first condition p\ > is not violated when 
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the condition A > is met. Hence only the violation of condition p 2 > gives rise to diffusion 
instability. Hence the condition for diffusive instability is given by 

H(k 2 ) = d x d 2 k A + \d 2 P*{B 1 - B 2 H*h{P\ D)) + d 1 H*(C 2 /P* - Ff 2 (H*, D 1 ))j k 2 + B<0. 

(4.7) 

H is quadratic in k 2 and the graph of H{k 2 ) = is a parabola. The minimum of H{k 2 ) occurs at 
k 2 = k 2 where 



d 2 P*{B 2 H*h{P\ D) - B t ) + d 1 H*(Ff 2 (H*, D x ) - C 2 /P* 



m 2d x d 2 L 

Consequently the condition for diffusive instability is H(k 2 n ). Therefore 
1 



Ad\d 2 



d 2 P*{B 2 H*h{P\ D) - B,) + d l H*(Ff 2 (H\ D x ) - C 2 /P* 



> 0. 



> B. 



(4.8) 



(4.9) 



Theorem 5. The criterion for diffusive instability for the model system is obtained by combining 
the result of Theorem 1, (14.81) and (14.91 ) and leading to the following condition: 

d 2 P*{B 2 H*f 1 {P*, D) - B t ) + d 1 H*(Ff 2 {H\ D 1 ) - C 2 /P* > 2{Bd 1 d 2 )l > 0. (4.10) 

In the following theorem, we are able to show the global stability behaviour of the positive 
equilibrium in the presence of diffusion. 



Theorem 6. (i) If the equilibrium E* of system (|3.1I) is globally asymptotically stable, the corre- 
sponding uniform steady state of system (|4.2I) is also globally asymptotically stable. 

( ii) If the equilibrium E* of system (13.11) is unstable even then the corresponding uniform steady 
state of system (14.21 ) can be made globally asymptotically stable by increasing the diffusion coeffi- 
cient di and d 2 to a sufficiently large value for strictly positive wave-number k > 0. 

Proof. For the sake of simplicity, let P(x, t) = P, H(x, t) = H . For x G [0, R] and t G [0, oo], let 
us define a functional: 



Vi 



R 



V(P, H)dx 



where V(P, H) is defined as Eq. (13 .6b . 

Differentiating V\ with respect to time t along the solutions of system (14.21) . we get 

" R (dVdP 
\dP~di 



dVi 
~~dT 



dV dH\ 

dH dt )' 



■ fdV dP 

' h — — )dx 



R dV , 
——dx 



o 



dt 



R -dV d 2 P dV d 2 H 

d\ TTZT — 7T + d 2 







OP dx 2 



' dH dx 2 



dx. 



Using the boundary condition (12.31) . we obtain 

dVi f R dV , , f R P*(P 2 + 3D 2 ) - 2D 2 P fdP\ 2 
——dx — di 



dt 



dt 



B 2 P 4 



R 



wH* /dH\* , 

-ht(&) dx - (4J1) 



From Eq. (14.1 II ). we note that if ^ < then in the interior of the positive quadrant of the 
PH— plane, and hence the first part of the theorem follows. We also note that if ^ < then ^ 
can be made negative by increasing d\ and d 2 a sufficiently large value, and hence the second part 
of the theorem follows. 
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4.2. Two dimensional case 

In two-dimensional case, the model system (12.11) can be written as 



B 1 P 2 - 



B 2 P 2 H 
P 2 +D 2 



a 2 p 

dx 2 



d 2 P 

W 1 



FH 2 



(4.12) 

- d 2 H , d 2 H \ 
H 2 +D 2 T" U 2 + ~dyr)- 

In this section, we show that the result of theorem 5 remain valid for two dimensional case. To 



prove this result we consider the following functional (Dubey and Hussain, 2000): 



V " = II V ^ H)dA 



where V(P, H) is defined as Eq. (13 .6b . 

Differentiating V 2 (t) with respect to time t along the solutions of system (14.121) , we obtain 



(W2 
dt 



where 



H dA ' 



Using Green's first identity in the plane 

J j FV 2 GdA = J F^-ds - J J (VF ■ VG)dA. 

q an n 

and under an analysis similar to Dubey and Hussain(2000), one can show that 



di 



d 9 



-di 



dV 
dH 



V 2 H dA 



d 2 V 
dP 2 

d 2 V 
dH 2 



,dP_ 

' dx 

' dx ' 



,dP, 
~dy ] 



dA < 0, 



dA < 0. 



(4.13) 



This shows that I 2 < 0. From the above analysis we note that if Ii < 0, then This implies that 
if in the absence of diffusion E* is globally asymptotically stable, then in the presence of diffusion 
E* will remain globally asymptotically stable. 

We also note that if ^ > 0, then I\ > 0. In such a case diffusion, E* will be unstable in the 
absence of diffusion. Even in this case by increasing di and d 2 to a sufficiently large value ^ can 
be made negative. This shows that if in the absence of diffusion E* is unstable, then in the presence 
of diffusion E* can be made stable by increasing diffusion coefficients to sufficiently large value. 
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5. Numerical simulations 

In this section, we perform numerical simulations to illustrate the results obtained in previous sec- 
tions. The dynamics of the model system (12.11) is studied with the help of numerical simulation, 
both in one and two dimensions, to investigate the spatiotemporal dynamics of the model sys- 
tem (I2.ll ). For the one-dimensional case, the plots (space vs. population densities) are obtained to 
study the spatial dynamics of the model system. The temporal dynamics is studied by observing 
the effect of time on space vs. density plot of prey populations. For the two-dimensional case, the 
spatial snapshots of prey densities are obtained at different time levels for different values of F 
and we have tried to study the spatiotemporal dynamics of the spatial model system. The rate of 
fish predation, F acting as forcing term or distributed control parameter in the model system (12.11 ) 
which is strictly nonnegative. Mathematically, we assume that F can be manipulated at every point 
in space and time. However, from practical point of view, we only have direct control of the net 
density of fish population at any instant. The parameters determining the shape of the type-Ill 
functional response of fish to zooplankton density are set more or less arbitrarily. They will in fact 
very much depend on the fish species, and also change significantly within a species with the size 
of the individuals. 

5.1. Spatiotemporal dynamics in One-dimensional system 

Now we start our insight into the spatiotemporal dynamics of the model system (12.11) in one- 
dimension with the following set of parameter values. The spatiotemporal dynamics of the system 
depends to a large extent on the choice of initial conditions. In a real aquatic ecosystem, the 
details of the initial spatial distribution of the species can be caused by spatially homogenous 
initial conditions. However, in this case, the distribution of species would stay homogenous for 
any time, and no spatial pattern can emerge. To get a nontrivial spatiotemporal dynamics, we 
have perturbed the homogenous distribution. We start with a hypothetical "constant-gradient" 
distribution (Malchow et al., 2008): 

P(x,0) = P*, H(x, 0) = H* + ex + 5, (5.1) 

where e and 5 are parameters. (P*, H*) = (1.8789, 1.3984) is the non-trivial equilibrium point of 
system (13.11 ), with fixed parameter set 

r = l t B 1 = 0.2, B 2 = 0.91, C x = 0.22, C 2 = 0.2, D 2 = 0.3, D x = 0.1, F = 0.1. 

It appears that the type of the system dynamics depends on e and 5. In Figs. 2 and 3, we 
show the population distribution over space at time t = 2000 (Fig. 2) and t = 1000, 2000(Fig.3), 
respectively. In Fig. 4, we illustrate the pattern formation about the above two cases, the system 
size is 3000, iteration numbers from 17800 to 19800, i.e., time t from 1780 to 1980. In Fig.2a 
and Fig.4a, s = —1.5 x 10~ 6 and 6 = 10^ 6 , the spatial distributions gradually vary in time, and 
the local temporal behavior of the dynamical variables P and H is strictly periodic following limit 
cycle of the non-spatial system, which is perhaps what is intuitively expected from system (14.11) . 
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there emerge a periodic wave pattern (Fig.4a, for P). However, in Fig. 2b and Fig. 4b, for a slightly 
different set of e and S, i.e., e = -1.5 x 1CT 3 and 5 = 10~ 6 , the dynamics of the system under- 
goes principal changes, and there emerge an oscillations wave pattern with a "flow" pattern around 
x = 1500 (c.f., Fig. 4b). In Fig. 3 and Fig.4c, e = —1.5 x 10 -2 and 5 = 10~ 3 , more different with 
the two cases above, the initial conditions (15.11) lead to the formation of strongly irregular "jagged" 
dynamics pattern inside a sub-domain of the system. There is a growing region of periodic travel- 
ling waves (x = 2000 to x = 2500) followed by oscillations with a large spatial wavelength. This 
is a typical scenario which precedes either the formation of spatiotemporal irregular oscillations or 
spatially homogeneous oscillations (c.f., Fig. 4c). 




I5CD 
Space 




(a) 



(b) 



Figure 2: Population distributions over space at t = 2000 obtained for the parameter values r = 
1,B 1 = 0.2, B 2 = 0.91, d = 0.22, C 2 = 0.2, D 2 = 0.3, A = 0.1, F = 0.1 and the initial 
conditions (ED with (a) e = -1.5 x 10" 6 and 5 = 10" 6 ; (b) e = -1.5 x 10~ 3 and 5 = 10~ 6 . 




500 1 000 1500 2000 2500 3000 500 1 000 1 500 2000 2500 3000 

Space Space 



Figure 3: Population distributions over space at (a) t = 1000 (b) t = 2000 obtained for the 
parameter values and the initial conditions (15.11 ) with £= -1.5 x 10~ 2 and 5 = 10~ 3 . The values 
of the other parameters are the same as given in Fig. 2. 
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(a) (b) (c) 

Figure 4: Spatiotemporal dynamical of system (14.11 ) with the initial conditions (15.11) . Parameters: 

(a) s = -1.5 x 10~ 6 , 5 = 10~ 6 ; (b) e = -1.5 x 1(T 3 , 5 = 1(T 6 ; (c) e = -1.5 x 10~ 2 , 5 = 1(T 3 . 

5.2. Spatiotemporal dynamics in Two-dimensional system 

In this subsection, we show the spatiotemporal dynamics in two-dimensional case of system (14.121) . 

5.2.1. Turing pattern 

First, we show Turing pattern of system (14.12t obtained by performing numerical simulations with 
initial and boundary conditions given in Eqs. (|2.2| )- (|2.31 ) with the parameter values given below in 
Eq. (27) with system size 200 x 200 for different values of F and time t. The values of the other 
parameters are fixed as 

r = l,Bi = 0.2, Bo = 0.91, d = 0.22, C 9 = 0.2, D 2 = 0.3, D x = 0.1, 

(5.2) 

d x = 0.005, d 2 = 1, Ax = Ay = 1/3, At = 1/40. 

Initially, the entire system is placed in the stationary state (P*, H*), and the propagation ve- 
locity of the initial perturbation is thus on the order of 5 x 10 -4 space units per time unit. And the 
system is then integrated for 10 5 or 2 x 10 5 time steps and some images saved. After the initial 
period during which the perturbation spreads, either the system goes into a time dependent state, or 
to an essentially steady state (time independent). Here, we show the distribution of phytoplankton 
P, for instance. 

Fig. 5 shows five typical Turing pattern of phytoplankton P in system (14.121) with parameters 
set (15.21) . From Fig. 5, one can see that values for the concentration P are represented in a color 
scale varying from blue (minimum) to red (maximum), and on increasing the control parameter 
F, the sequences spots — > spot-stripe mixtures— s> stripes— > hole-stripe mixtures holes pattern is 
observed. When F > 0.077, other parameters are the same as above case, wave pattern emerges 
and are presented in Fig. 6. 

5.2.2. Spiral wave pattern 

Thanks to the insightful work of Medvinsky et al. (2002), we have studied the spiral wave pattern 
for two different set of initial condition discussed in Eqs. ( 15.31) and ( 15.41) respectively. In these two 
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Figure 5: Typical Turing patterns of P in the system (14.121) with fixed parameters set (I5.2I ). (a) 
spots pattern, F = 0.0001; (b) spot-stripe mixtures pattern, F = 0.01; (c) stripes pattern, F = 
0.035; (d) hole-stripe mixtures, F = 0.065; (e) holes pattern, F = 0.0765. Iterations: Pattern-(a) 
and(e): 2 x 10 5 , others: 10 5 . 
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Figure 6: Wave pattern obtained with system (14.121 ) with F = 0.085. Iterations: (a) 0, (b) 10000, 
(c) 25000, (d) 40000. 
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cases, we employ Ax = Ay = 1, At = 1/3 and the system size is 900 x 300. Other parameters 
set is fixed as (15.2b . 
The first case is: 

P(x, y, 0) = P* - ex(x - 0.2y - 225)(x - 0.2y - 675), 

H(x, y, 0)=H*- e 2 {x - 450) - e 3 {y - 150). (5 ' 3) 

where e 1 = 2 x 10~ 7 , e 2 = 3 x 10~ 5 , e 3 = 1.2 x 10~ 4 . 

The initial conditions are deliberately chosen to be unsymmetrical in order to make any in- 
fluence of the comers of the domain more visible. Snapshots of the spatial distribution arising 
from (Q are shown in Fig. 7 for t = 0, 2500, 5000, 10000. Fig.7a shows that for the system (I4T2T) 
with initial conditions (15.31 ), the formation of the irregular patchy structure can be preceded by the 
evolution of a regular spiral spatial pattern. Note that the appearance of the spirals is not induced 
by the initial conditions. The center of each spiral is situated in a critical point (x cr , y cr ), where 
U(x cr ,y cr ) = P*, V(x cr ,y cr ) = H*. The distribution (15.31) contains two such points; for other 
initial conditions, the number of spirals may be different. After the spirals form (Fig. 7b), they 
grow slightly for a certain time, their spatial structure becoming more distinct (Figs. 7c and 7d). 




Figure 7: Spiral pattern of P in system (14.121 ) with special initial condition (15.31) and F = 0.1. 
Time: (a) 0, (b) 2500, (c) 5000, (d) 10000. 

In the second case, the initial conditions describe a phytoplankton (prey) patch placed into a 
domain with a constant-gradient zooplankton (predator) distribution: 

P(x, y, 0) = P* - e 1 {x - 180)(x - 720) - s 2 (y - 90) (y - 210), 

(5.4) 

H(x, y, 0) = H* - e 3 {x - 450) - e A (y - 135). 

where £i = 2x 10~ 7 , e 2 = 6 x 10~ 7 , e 3 = 3 x 10~ 5 , e 4 = 6 x 10" 5 . 

Fig. 8 shows the snapshots of the spatial distribution arising from (15.31 ) for t — 0, 1000, 4000, 9000. 
Although the dynamics of the system preceding the formation of patchy spatial structure is some- 
what less regular than in the previous case, it follows a similar scenario. Again the spirals appear, 
with their centers located in the vicinity of critical points (Figs. 8c and 8d). 



15 



R. K. Upadhyay et al. 



Spatiotemporal dynamics in a spatial plankton system 




3.5 



2.5 



1,5 



(a) 



(b) 




(c) 



(d) 



Figure 8: Spiral pattern of P in system (I4.12I ) with special initial condition (I5.4I) and F = 0.1. 
Time: (a) 0, (b) 1000, (c) 4000, (d) 9000. 

6. Discussions and Conclusions 

In this paper, we have made an attempt to provide a unified framework to understand the com- 
plex dynamical patterns observed in a spatial plankton model for phytoplankton-zooplankton-fish 
interaction with Holling type III functional response. Since most marine mammals and boreal 
fish species are considered to be generalist predators (including feeders, grazers, etc.), a sigmoidal 
(type III) functional response might be more appropriate (Magnusson and Palsson, 1991). This 
has been proposed as being likely for minke whales, harp seals, and cod in the Barents Sea where 
switching between prey species has been hypothesized (Nilssen et al., 2000; Schweder et al., 2000). 
Sigmoidal functional response arising due to heterogeneity of space is less addressed in the litera- 
ture. It is well known that real vertical distribution of plankton is highly heterogeneous, especially 
during plankton blooms and this can seriously affect the behavior of plankton models (Poggiale et 
al. 2008). Morozov (2010) shows that emergence of Holling type III in plankton system is due to 
mechanisms different from those well known in the literature (e.g. food search learning by preda- 
tor, existence of alternative food, refuse for prey). It has certainly been useful in both theoretical 
investigations and practical applications (Turchin, 2003). Movement of phytoplankton and zoo- 
plankton population with different velocities can give rise to spatial patterns (Malchow, 2000). We 
have studied the reaction-diffusion model in both one and two dimensions and investigated their 
stability conditions. By analyzing the model system (12.11) . we observed that two points of equi- 
librium exist for the model system. The non-trivial equilibrium state of prey-predator coexistence 
is locally as well as globally asymptotically stable under a fixed region of attraction when certain 
conditions are satisfied. It has been shown that the instability observed in the system is diffusion- 
driven under the condition (I4.101 ). The expression for critical wave-number has been obtained. It 
has been observed that the solution of the model system converges faster to its equilibrium in case 
of diffusion in two-dimensional space as compared to the diffusion in one dimensional space. 
In the numerical simulation, we adopt the rate of fish predation F in the system (12.11) as the con- 
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trol parameter which is strictly nonnegative. In the two-dimensional case, on increasing the value 
of F, the sequences spots— i>spot-stripe mixtures— ^stripes— > hole-stripe mixtures— wholes— >wave 
pattern is observed. For the sake of learning the wave pattern further, we show the time evolution 
process of the pattern formation with two special initial conditions and find spiral wave pattern 
emerge. That is to say, our two-dimensional spatial patterns may indicate the vital role of phase 
transience regimes in the spatiotemporal organization of the phytoplankton and zooplankton in the 
aquatic systems. It is also important to distinguish between "intrinsic" patterns, i.e., patterns aris- 
ing due to trophic interaction such as those considered above, and "forced" patterns induced by 
the inhomogeneity of the environment. The physical nature of the environmental heterogeneity, 
and thus the value of the dispersion of varying quantities and typical times and lengths, can be 
essentially different in different cases. 
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Appendix: Proof of Theorem 1 



The Jacobian matrix of model (13.11) corresponding to the positive equilibrium E*(P*, H*) is given 
by 

fen 612 
hi hi 



J 



where 



hi=r- 1B V P* 
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Since E*(P*, H*) is a solution of Eq. (13.21) . From the first equation of Eq. (13 .2b . we obtain 
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From the second equation of Eq. (I3.21 ). we have 
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where f 2 (H*,D 1 )= ffi^l 

Now the characteristic equation corresponding to the matrix J is given by 

A 2 - (hx + 6 22 )A + (611622 - 612621) or A 2 + A\ + B = 0, 



where A and B are defined in (13.31) . By Routh-Hurwitz criteria, all eigenvalues of J will have 
negative real parts if and only if A > 0, B > 0. Thus, the theorem follows. 
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